#include "NeutrinoTools.h"

TString getTitleFromVar(const TString var)
{
  TString tit(var);
  tit.ReplaceAll("fMomErr/fMuonRec->P()","#sigma_{p}/p rec");
  tit.ReplaceAll("fNeutrinoSim->E()","E_{#nu sim} (GeV)");
  tit.ReplaceAll("fMuonSim->P()","p_{#mu sim} (GeV/c)");
  tit.ReplaceAll("fMuonSim->Pt()","p_{t}^{#mu sim int} (GeV/c)");
  tit.ReplaceAll("fMuonRec->Pt()","p_{t}^{#mu rec int} (GeV/c)");
  //tit.ReplaceAll("fMuonPtGuess->Mag()","p_{t}^{#mu rec, #nu sim} (GeV/c)");
  tit.ReplaceAll("fMuonSimPt->Mag()","p_{t}^{#mu sim} (GeV/c)");
  tit.ReplaceAll("fMuonRec->P()","p_{#mu rec} (GeV/c)");
  tit.ReplaceAll("fMuonRecPt->Mag()","p_{t}^{#mu rec} (GeV/c)");
  tit.ReplaceAll("fQ2Sim","Q^{2}_{sim} (GeV^{2}/c^{4})");
  tit.ReplaceAll("fMultiplicity","N_{all rec}");
  tit.ReplaceAll("fMuonRecFlightPath->Pt()","lever arm (m)");
  tit.ReplaceAll("fMuonRecFlightPath->Mag()","straight path length (m)");
  tit.ReplaceAll("(fMuonCharge[1]==-1)*2+(fMuonCharge[0]==-1)","TPC charge (11,-11,1-1,-1-1)");
  tit.ReplaceAll("fMuonRecNhits","selmu_tpc_nhits");
  tit.ReplaceAll("fMuonTypeSim","1:e, 2:#mu, 3:K_{L}, 4:#pi^{#pm}, 5:K^{#pm}, 6:p");
  tit.ReplaceAll("fNeutrinoParentType","1:e, 2:#mu, 3:K_{L}, 4:#pi^{#pm}, 5:K^{#pm}, 6:p");
  tit.ReplaceAll("fMuonAlphaRec","#theta_{#mu rec} (deg)");
  tit.ReplaceAll("fMuonAlphaSim","#theta_{#mu sim} (deg)");
  tit.ReplaceAll("fNeutrinoRec->Theta()*TMath::RadToDeg()","#theta_{#nu rec} (deg)");
  tit.ReplaceAll("fNeutrinoSim->Theta()*TMath::RadToDeg()","#theta_{#nu sim} (deg)");
  tit.ReplaceAll("fMuonRecEndPos->X()","endposX (m)");
  tit.ReplaceAll("fMuonRecEndPos->Y()","endposY (m)");
  tit.ReplaceAll("fMuonRecEndPos->Z()","endposZ (m)");
  tit.ReplaceAll("fMuonRecVertex->X()","vertexX (m)");
  tit.ReplaceAll("fMuonRecVertex->Y()","vertexY (m)");
  tit.ReplaceAll("fMuonRecVertex->Z()","vertexZ (m)");
  tit.ReplaceAll("fMuonRec->Theta()*TMath::RadToDeg()","#theta_{#mu rec lab} (deg)");
  tit.ReplaceAll("(fNeutMode==1||fNeutMode==2)","CCQE");
  tit.ReplaceAll("fMuonRecChi2/fMuonRecNDOF","#chi^{2}_{track}/NDOF");
}

plotMom(){
  TString tag="mom";
  gStyle->SetOptStat(0);

  t0=new TH1D("t0","",100,0.1,20); NeutrinoTools::BinLog(t0->GetXaxis()); t0->SetTitle("red: no cut, blue: M>1, green: M>1 & mu_P > 1 GeV/c ;Q^{2}_{sim} (GeV^2);Counts"); t0->SetLineColor(kRed); t0->Sumw2(); t0->SetLineWidth(2);
  t1=new TH1D("t1","",100,0.1,20); NeutrinoTools::BinLog(t1->GetXaxis()); t1->SetLineColor(kBlue); t1->Sumw2();t1->SetLineWidth(2);
  t2=new TH1D("t2","",100,0.1,20); NeutrinoTools::BinLog(t2->GetXaxis()); t2->SetLineColor(kGreen); t2->Sumw2();t2->SetLineWidth(2);
  
  tree->Draw("fQ2Sim>>t0");
  tree->Draw("fQ2Sim>>t1","fMultiplicity>1");
  tree->Draw("fQ2Sim>>t2","fMultiplicity>1 && fMuonRec.P()>1");

  e10=(TH1D*)t1->Clone("e10"); e10->Divide(t1,t0,1,1,"B"); e10->SetTitle("blue: M>1, green: M>1 & mu_P > 1 GeV/c, black: green/blue; Q^{2}_{sim} (GeV^2); #epsilon"); 
  e20=(TH1D*)t2->Clone("e20"); e20->Divide(t2,t0,1,1,"B");
  e21=(TH1D*)t2->Clone("e21"); e21->Divide(t2,t1,1,1,"B"); e21->SetLineColor(kBlack);
  
  c1=new TCanvas;
  gPad->SetLogx();
  gPad->SetLogy();
  t0->Draw();
  t1->Draw("same");
  t2->Draw("same");
  c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_count.png");
  c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_count.eps");

  c2=new TCanvas;
  gPad->SetLogx();
  e10->Draw("e");
  e20->Draw("same e");
  e21->Draw("same e");
  c2->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_eff.png");
  c2->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_eff.eps");
}

plotPt(){
  TString tag="pt";
  gStyle->SetOptStat(0);

  t0=new TH1D("t0","",100,0.1,20); NeutrinoTools::BinLog(t0->GetXaxis()); t0->SetTitle("red: no cut, blue: M>1, green: M>1 & mu_pT > 1 GeV/c ;Q^{2}_{sim} (GeV^2);Counts"); t0->SetLineColor(kRed); t0->Sumw2();t0->SetLineWidth(2);
  t1=new TH1D("t1","",100,0.1,20); NeutrinoTools::BinLog(t1->GetXaxis()); t1->SetLineColor(kBlue); t1->Sumw2();t1->SetLineWidth(2);
  t2=new TH1D("t2","",100,0.1,20); NeutrinoTools::BinLog(t2->GetXaxis()); t2->SetLineColor(kGreen); t2->Sumw2();t2->SetLineWidth(2);

  tree->Draw("fQ2Sim>>t0");
  tree->Draw("fQ2Sim>>t1","fMultiplicity>1");
  tree->Draw("fQ2Sim>>t2","fMultiplicity>1 && fMuonRecPt.Mag()>1");

  e10=(TH1D*)t1->Clone("e10"); e10->Divide(t1,t0,1,1,"B"); e10->SetTitle("blue: M>1, green: M>1 & mu_pT > 1 GeV/c, black: green/blue; Q^{2}_{sim} (GeV^2); #epsilon");
  e20=(TH1D*)t2->Clone("e20"); e20->Divide(t2,t0,1,1,"B");
  e21=(TH1D*)t2->Clone("e21"); e21->Divide(t2,t1,1,1,"B"); e21->SetLineColor(kBlack);

  c1=new TCanvas;
  gPad->SetLogx();
  gPad->SetLogy();
  t0->Draw();
  t1->Draw("same");
  t2->Draw("same");
  c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_count.png");
  c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_count.eps");

  c2=new TCanvas;
  gPad->SetLogx();
  e10->Draw("e");
  e20->Draw("same e");
  e21->Draw("same e");
  c2->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_eff.png");
  c2->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_eff.eps");
}

plotPtGuess(){
  TString tag="PtGuess";
  gStyle->SetOptStat(0);

  t0=new TH1D("t0","",100,0.1,20); NeutrinoTools::BinLog(t0->GetXaxis()); t0->SetTitle("red: no cut, blue: M>1, green: M>1 & mu_PtGuess > 1 GeV/c ;Q^{2}_{sim} (GeV^2);Counts"); t0->SetLineColor(kRed); t0->Sumw2();t0->SetLineWidth(2);
  t1=new TH1D("t1","",100,0.1,20); NeutrinoTools::BinLog(t1->GetXaxis()); t1->SetLineColor(kBlue); t1->Sumw2();t1->SetLineWidth(2);
  t2=new TH1D("t2","",100,0.1,20); NeutrinoTools::BinLog(t2->GetXaxis()); t2->SetLineColor(kGreen); t2->Sumw2();t2->SetLineWidth(2);

  tree->Draw("fQ2Sim>>t0");
  tree->Draw("fQ2Sim>>t1","fMultiplicity>1");
  tree->Draw("fQ2Sim>>t2","fMultiplicity>1 && fMuonPtGuess>1");

  e10=(TH1D*)t1->Clone("e10"); e10->Divide(t1,t0,1,1,"B"); e10->SetTitle("blue: M>1, green: M>1 & mu_PtGuess > 1 GeV/c, black: green/blue; Q^{2}_{sim} (GeV^2); #epsilon");
  e20=(TH1D*)t2->Clone("e20"); e20->Divide(t2,t0,1,1,"B");
  e21=(TH1D*)t2->Clone("e21"); e21->Divide(t2,t1,1,1,"B"); e21->SetLineColor(kBlack);

  c1=new TCanvas;
  gPad->SetLogx();
  gPad->SetLogy();
  t0->Draw();
  t1->Draw("same");
  t2->Draw("same");
  c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_count.png");
  c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_count.eps");

  c2=new TCanvas;
  gPad->SetLogx();
  e10->Draw("e");
  e20->Draw("same e");
  e21->Draw("same e");
  c2->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_eff.png");
  c2->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_eff.eps");
}


plotPtInv(){
  TString tag="ptinv";

  gStyle->SetOptStat(0);

  t0=new TH1D("t0","",100,0.1,20); NeutrinoTools::BinLog(t0->GetXaxis()); t0->SetTitle("red: no cut, blue: mu_pT>1, green: M>1 & mu_pT > 1 GeV/c ;Q^{2}_{sim} (GeV^2);Counts"); t0->SetLineColor(kRed); t0->Sumw2();t0->SetLineWidth(2);
  t1=new TH1D("t1","",100,0.1,20); NeutrinoTools::BinLog(t1->GetXaxis()); t1->SetLineColor(kBlue); t1->Sumw2();t1->SetLineWidth(2);
  t2=new TH1D("t2","",100,0.1,20); NeutrinoTools::BinLog(t2->GetXaxis()); t2->SetLineColor(kGreen); t2->Sumw2();t2->SetLineWidth(2);

  tree->Draw("fQ2Sim>>t0");
  tree->Draw("fQ2Sim>>t1","fMuonRecPt.Mag()>1");
  tree->Draw("fQ2Sim>>t2","fMultiplicity>1 && fMuonRecPt.Mag()>1");

  e10=(TH1D*)t1->Clone("e10"); e10->Divide(t1,t0,1,1,"B"); e10->SetTitle("blue: mu_pT>1, green: M>1 & mu_pT > 1 GeV/c, black: green/blue; Q^{2}_{sim} (GeV^2); #epsilon");
  e20=(TH1D*)t2->Clone("e20"); e20->Divide(t2,t0,1,1,"B");
  e21=(TH1D*)t2->Clone("e21"); e21->Divide(t2,t1,1,1,"B"); e21->SetLineColor(kBlack);

  c1=new TCanvas;
  gPad->SetLogx();
  gPad->SetLogy();
  t0->Draw();
  t1->Draw("same");
  t2->Draw("same");
  c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_count.png");
  c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_count.eps");

  c2=new TCanvas;
  gPad->SetLogx();
  e10->Draw("e");
  e20->Draw("same e");
  e21->Draw("same e");
  c2->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_eff.png");
  c2->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_eff.eps");
}

plotPtGuessInv()
{
  TString tag="PtGuessinv";

  gStyle->SetOptStat(0);

  t0=new TH1D("t0","",100,0.1,20); NeutrinoTools::BinLog(t0->GetXaxis()); t0->SetTitle("red: no cut, blue: mu_PtGuess>1, green: M>1 & mu_PtGuess > 1 GeV/c ;Q^{2}_{sim} (GeV^2);Counts"); t0->SetLineColor(kRed); t0->Sumw2();t0->SetLineWidth(2);
  t1=new TH1D("t1","",100,0.1,20); NeutrinoTools::BinLog(t1->GetXaxis()); t1->SetLineColor(kBlue); t1->Sumw2();t1->SetLineWidth(2);
  t2=new TH1D("t2","",100,0.1,20); NeutrinoTools::BinLog(t2->GetXaxis()); t2->SetLineColor(kGreen); t2->Sumw2();t2->SetLineWidth(2);

  tree->Draw("fQ2Sim>>t0");
  tree->Draw("fQ2Sim>>t1","fMuonPtGuess>1");
  tree->Draw("fQ2Sim>>t2","fMultiplicity>1 && fMuonPtGuess>1");

  e10=(TH1D*)t1->Clone("e10"); e10->Divide(t1,t0,1,1,"B"); e10->SetTitle("blue: mu_PtGuess>1, green: M>1 & mu_PtGuess > 1 GeV/c, black: green/blue; Q^{2}_{sim} (GeV^2); #epsilon");
  e20=(TH1D*)t2->Clone("e20"); e20->Divide(t2,t0,1,1,"B");
  e21=(TH1D*)t2->Clone("e21"); e21->Divide(t2,t1,1,1,"B"); e21->SetLineColor(kBlack);

  c1=new TCanvas;
  gPad->SetLogx();
  gPad->SetLogy();
  t0->Draw();
  t1->Draw("same");
  t2->Draw("same");
  c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_count.png");
  c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_count.eps");

  c2=new TCanvas;
  gPad->SetLogx();
  e10->Draw("e");
  e20->Draw("same e");
  e21->Draw("same e");
  c2->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_eff.png");
  c2->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"_eff.eps");
}

/*

.L plot.C
double tol=0.05; TString tag = "muonRecPtNeutrinoE"; Double_t x0 = 0.09, x1 = 11; Bool_t klogx=1; while(tol<=0.21){ plotFlux(tol, tag, x0, x1, klogx); tol+=0.05;}
double tol=0.05; TString tag = "AlphaRecNeutrinoE"; Double_t x0 = 0, x1 = 180; Bool_t klogx=0; while(tol<=0.21){ plotFlux(tol, tag, x0, x1, klogx); tol+=0.05;}
double tol=0.05; TString tag = "muonPRecNeutrinoE"; Double_t x0 = 0.09, x1 = 15; Bool_t klogx=1; while(tol<=0.21){ plotFlux(tol, tag, x0, x1, klogx); tol+=0.05;}
double tol=0.05; TString tag = "MultiplicityNeutrinoE"; Double_t x0 = 0, x1 = 10; Bool_t klogx=0; while(tol<=0.21){ plotFlux(tol, tag, x0, x1, klogx); tol+=0.05;}

*/

plotFlux(const Double_t tol = 0.05, const TString tag="muonRecPtNeutrinoE", const Double_t x0=0.09, const Double_t x1=11, const Bool_t klogx=1)
{
  gStyle->SetOptStat(0);

  hh=(TH2D*) gDirectory->Get(tag);
  if(!hh){
    printf("no hist %s\n", tag.Data()); return;
  }

  TH2D *hmap=0x0;
  TH2D *hdiff=0x0;

  TGraphAsymmErrors *gr = NeutrinoTools::GetFluxMap(tol, hh, hmap, hdiff);
  
  c1=new TCanvas; gPad->SetLogx(); gPad->SetLogy(klogx);  gPad->SetLogz(); hh->Draw("colz"); c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"original.png");
  c1=new TCanvas; gPad->SetLogx(); gPad->SetLogy(klogx);  hmap->SetMinimum(0); hmap->SetMaximum(2); hmap->Draw("colz"); c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"map1-2.png");
  c1=new TCanvas; gPad->SetLogx(); gPad->SetLogy(klogx);  hmap->SetMinimum(0.5); hmap->SetMaximum(1.5); hmap->Draw("colz"); c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"map0.5-1.5.png");
  
  c1=new TCanvas; gPad->SetLogx(); gPad->SetLogy(klogx);  hdiff->SetMinimum(0); hdiff->SetMaximum(1); hdiff->Draw("colz"); c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"diff0-1.png");
  c1=new TCanvas; gPad->SetLogx(); gPad->SetLogy(klogx);  hdiff->SetMinimum(0); hdiff->SetMaximum(0.1); hdiff->Draw("colz"); c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"diff0-0.1.png");
  
  c1=new TCanvas; gPad->SetGrid();  gPad->SetLogx(klogx); gPad->SetLogy(); gr->SetMarkerSize(1); gr->SetMarkerStyle(24); gr->GetXaxis()->SetLimits(x0, x1); gr->GetYaxis()->SetRangeUser(0.2, 21); gr->Draw("alp"); c1->Print(Form("outplot/%scurveTol%.2f.png", tag.Data(), tol));

}

plotEMode(const Bool_t khigh=kTRUE)
{
  TCanvas *c1=new TCanvas;

  gStyle->SetOptStat(0);

  const Int_t mm = 30;
  const TString cut0 = "fMuonNSim>0 &&";
  const TString ecut = khigh? "fNeutrinoSim->E()>5.318&&fNeutrinoSim->E()<11.774" : "fNeutrinoSim->E()>0.490&&fNeutrinoSim->E()<1.085";

  TH1D * hh1 = new TH1D("hh1","",mm,-0.5, mm-0.5); hh1->SetLineWidth(4); 
  TH1D * hh2 = new TH1D("hh2","",mm,-0.5, mm-0.5); hh2->SetLineWidth(4); 

  //mode<0 for anti-neutrino
  tree->Draw("abs(fNeutMode)>>hh1",cut0+ecut+" && fQ2Sim<0.9"); hh1->SetLineColor(kRed);  hh1->SetTitle(ecut+", red: Q^{2}<0.9 GeV^{2}, blue: Q^{2}>1 GeV^{2}; Neut Mode; Counts");  
  tree->Draw("abs(fNeutMode)>>hh2",cut0+ecut+" && fQ2Sim>1");   hh2->SetLineColor(kBlue); 

  printf("hh1 ou %.0f hh2 ou %.0f\n", hh1->GetBinContent(0)+hh1->GetBinContent(hh1->GetNbinsX()), hh2->GetBinContent(0)+hh2->GetBinContent(hh2->GetNbinsX()));

  const Double_t hm = khigh? hh2->GetMaximum() : hh1->GetMaximum();

  hh1->SetMaximum( hm*1.1); hh1->SetMinimum(0);hh1->Draw();hh2->Draw("same");
  c1->Print(Form("outplot/neutmodelogy%dkhigh%d.png", gPad->GetLogy(), khigh));
  c1->Print(Form("outplot/neutmodelogy%dkhigh%d.eps", gPad->GetLogy(), khigh));
  
  gPad->SetLogy(1);
  hh1->SetMaximum( hm*2); hh1->SetMinimum(0.5);hh1->Draw();hh2->Draw("same");
  c1->Print(Form("outplot/neutmodelogy%dkhigh%d.png", gPad->GetLogy(), khigh));
  c1->Print(Form("outplot/neutmodelogy%dkhigh%d.eps", gPad->GetLogy(), khigh));
}

plotptMode()
{
  TCanvas *c1=new TCanvas;

  gStyle->SetOptStat(0);

  const Int_t mm = 30;
  TH1D * hh1 = new TH1D("hh1","",mm,-0.5, mm-0.5); hh1->SetLineWidth(5);
  TH1D * hh2 = new TH1D("hh2","",mm,-0.5, mm-0.5); hh2->SetLineWidth(5);
  TH1D * hh3 = new TH1D("hh3","",mm,-0.5, mm-0.5); hh3->SetLineWidth(4);
  TH1D * hh4 = new TH1D("hh4","",mm,-0.5, mm-0.5); hh4->SetLineWidth(3);

  tree->Draw("abs(fNeutMode)>>hh1","fMuonRecPt->Mag()<1"); hh1->SetLineColor(kRed);  
  tree->Draw("abs(fNeutMode)>>hh2","fMuonRecPt->Mag()>1"); hh2->SetLineColor(kBlue); 
  tree->Draw("abs(fNeutMode)>>hh3","fMuonRecPt->Mag()>1.5"); hh3->SetLineColor(kGreen); 
  tree->Draw("abs(fNeutMode)>>hh4","fMuonRecPt->Mag()>2"); hh4->SetLineColor(kMagenta); 
  const Double_t sf3 = hh2->GetBinContent(27)/hh3->GetBinContent(27);
  hh3->Scale(sf3);

  const Double_t sf4 = hh2->GetBinContent(27)/hh4->GetBinContent(27);
  hh4->Scale(sf4);

  const TString tit=Form("red: p_{t}^{#mu, rec}<1 GeV/c, blue: p_{t}^{#mu, rec}>1 GeV/c, green: #times %.1f p_{t}^{#mu, rec}>1.5 GeV/c, magenta: #times %.1f p_{t}^{#mu, rec}>2 GeV/c; Neut Mode; Counts", sf3, sf4);  
  hh1->SetTitle(tit);hh2->SetTitle(tit);

  printf("OU: hh1 %.0f hh2 %.0f hh3 %.0f hh4 %0.f\n", hh1->GetBinContent(0)+hh1->GetBinContent(hh1->GetNbinsX()), hh2->GetBinContent(0)+hh2->GetBinContent(hh2->GetNbinsX()), hh3->GetBinContent(0)+hh3->GetBinContent(hh3->GetNbinsX()), hh4->GetBinContent(0)+hh4->GetBinContent(hh4->GetNbinsX()));

  const Double_t hm = hh1->GetMaximum();

  hh2->Draw(); hh3->Draw("same"); hh4->Draw("same");
  c1->Print(Form("outplot/PtModelogy%dk.png", gPad->GetLogy()));
  c1->Print(Form("outplot/PtModelogy%dk.eps", gPad->GetLogy()));

  gPad->SetLogy(1);
  hh1->SetMaximum( hm*2); hh1->SetMinimum(0.5);hh1->Draw();hh2->Draw("same"); hh3->Draw("same"); hh4->Draw("same");
  c1->Print(Form("outplot/PtModelogy%dk.png", gPad->GetLogy()));
  c1->Print(Form("outplot/PtModelogy%dk.eps", gPad->GetLogy()));


}


plotNuE(const TString tag)
{
  ll=lout;

  TDirectory *tmpdir = new TDirectory("tmpdir","tmpdir");
  for(Int_t ii=0; ii<ll->GetEntries(); ii++){
    htmp = ll->At(ii);
    //printf("loading %d: %s\n", ii, htmp->GetName());
    tmpdir->Add(htmp);
  }

  tmpdir->cd();

  c1=new TCanvas;
  gStyle->SetOptStat(0);
  NeutrinoE->SetTitle((tag.Contains("/")?tag:tag.Prepend("outplot/"))+";E_{#nu, sim}; Counts");
  gPad->SetLogx();
  NeutrinoE->Draw();

  c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"nuE.png");
  c1->Print((tag.Contains("/")?tag:tag.Prepend("outplot/"))+"nuE.eps");

  delete tmpdir;
}

plotResPtGuess()
{
  gROOT->LoadMacro("style.cxx+");

  style::SetGlobalStyle();

  style::fgkYTitleOffset = 1.4;
  style::fgkXTitleOffset = 1.4;

  TCanvas *c1= new TCanvas;
  style::PadSetup(c1);

  Double_t lm=0.15;
  Double_t rm=0.13;
  Double_t tm=0.03;
  Double_t bm=0.17;

  c1->SetLeftMargin(lm);
  c1->SetRightMargin(rm);
  c1->SetTopMargin(tm);
  c1->SetBottomMargin(bm);

  TH2D * hh = new TH2D ("hh","",50,0.1, 2.5,50,-1,1); NeutrinoTools::BinLog(hh->GetXaxis());
  hh->SetTitle(";p_{t}^{#mu sim} (GeV/c);p_{t}^{#mu rec}/p_{t}^{#mu sim}-1");
  style::ResetStyle(hh);
  tree->Draw("fMuonPtGuess->Mag()/fMuonSimPt->Mag()-1:fMuonSimPt->Mag()>>hh","fMuonNSim>0 && fMuonNRec>0 &&  fMultiplicity>=3");
  TH2D *hnor = NeutrinoTools::NormalHist(hh, 5, 1);
  style::ResetStyle(hnor);
  gPad->SetLogx();

  hh->Draw("colz");
  //c1->Print("outplot/PtGuessRes.eps");
  c1->Print("outplot/PtGuessRes.png");

  hnor->Draw("colz");
  //c1->Print("outplot/PtGuessNorRes.eps");
  c1->Print("outplot/PtGuessNorRes.png");
}

//====================================================================================
//====================================================================================

TString GetPIDCut()
{
  //return "( fMuonRecNhits[0]<80 && fMuonRecFlightPath->Mag()>2.5 && fMuonCharge[0]==-1 && fMuonCharge[1]==-1 && fMultiplicity<4)";
  return "( fMuonRecNhits[0]<80 && fMuonRecFlightPath->Mag()>3.0 )";
}

TString GetPtResCut()
{
  //don't cut on p, at least not on lower end, it is physics; too large P is mostly fluctuation
  //return "( fMuonRec->P()<3.5 && fMuonRecNhits[0]>60 && fMuonRecFlightPath->Pt()>1.5  )";
  return "( 0.25<(fMuonRecChi2/fMuonRecNDOF) && (fMuonRecChi2/fMuonRecNDOF)<4 && fMuonRecNhits[0]>60 && fMuonRecFlightPath->Pt()>1.5 )";
}


TString GetQACut()
{
  return Form("( %s && %s )", GetPIDCut().Data(), GetPtResCut().Data());
}

TString GetModeCut(const Int_t imode, TString & name)
{
  const TString modecuts[]={"1",
                            "fNeutrinoType==14 && ( fNeutMode==1 || fNeutMode==2 )",
                            "fNeutrinoType==14 && ( fNeutMode==11 || fNeutMode==12 || fNeutMode==13 || fNeutMode== 16 || fNeutMode==17 || fNeutMode== 22 || fNeutMode== 23 )",
                            "fNeutrinoType==14 && ( fNeutMode==21 )", 
                            "fNeutrinoType==14 && ( fNeutMode==26 )",
                            //"fNeutrinoType!=14 && ( abs(fNeutMode)<30 )"
                            "(fNeutrinoType!=14 || abs(fNeutMode)>30)"
  };
  const TString modenames[]={"incl","#nu_{#mu}CCQE","#nu_{#mu}CC#Delta","#nu_{#mu}CCm#pi","#nu_{#mu}CCDIS","#slash{#nu_{#mu}} or NC"};

  name=modenames[imode];
  return Form("(%s)",modecuts[imode].Data());
}

void plotStack(const TString fin, const TString xvar, const TString cvar, Bool_t kqacut0, Bool_t kosci0=kTRUE)
{
  const TString pwd=gSystem->pwd();

  if(pwd.Contains("SK")){
    if(cvar=="PID"||cvar.Contains("Parent") || cvar.Contains("Nhits") || xvar.Contains("Nhits") || cvar.Contains("Flight") || xvar.Contains("Flight") || cvar.Contains("Charge") || xvar.Contains("Charge")  || cvar.Contains("LeptonType") || xvar.Contains("LeptonType")  ){
      return;
    }
  }

  if(pwd.Contains("SK")){
    kqacut0=kFALSE;
  }
  else{
    kosci0=kFALSE;
  }

  if(pwd.Contains("ND280")){
    if(!(cvar=="CCpurity" && (xvar=="fNeutMode" || xvar=="fNeutrinoType"))){
      if((cvar=="CCeff" || cvar=="CCpurity") && kqacut0==0){
        return;
      }
    }
  }

  if(cvar=="Q2PPt" && xvar=="fQ2Sim"){
    return;
  }

  if(cvar=="EnuPPt" && xvar=="fNeutrinoSim->E()"){
    return;
  }

//  if(cvar=="CCpurity" && (xvar=="fNeutMode" || xvar=="fNeutrinoType")){
//    return;
//  }

  if(cvar=="CCeff" && (xvar=="fNeutrinoType")){
    return;
  }

  if(cvar=="PID" && xvar=="fMuonTypeSim"){
    return;
  }

  const Bool_t kqacut = kqacut0;
  const Bool_t kosci = kosci0;

  gROOT->LoadMacro("style.cxx+");

  TFile::Open(fin);

  TCanvas *c1=new TCanvas;
  style::SetGlobalStyle();
  style::PadSetup(c1);

  Double_t lm=0.15;
  Double_t rm=0.03;
  Double_t tm=0.08;
  Double_t bm=0.17;

  c1->SetLeftMargin(lm);
  c1->SetRightMargin(rm);
  c1->SetTopMargin(tm);
  c1->SetBottomMargin(bm);

  style::fgkYTitleOffset = 1.4;
  style::fgkXTitleOffset = 1.4;

  gStyle->SetOptStat(0);

  gPad->SetLogx();

  //-----------------------------------------------------------

  Int_t xnbin = 50;
  Double_t xmin = 1e-3;
  const TString pwd=gSystem->pwd();
  Double_t ptmax = 3;
  if(pwd.Contains("ND280") && xvar.Contains("Rec")){
    ptmax = 3;//kqacut?3:20;
  }
  Double_t xmax = xvar.Contains("Pt")?ptmax:20;


  const TString cutCCRec = Form("(fMuonNRec > 0 %s)", kqacut? Form(" && %s",(GetQACut()).Data()): "");
  const TString conCCRec = Form("#mu^{-} rec %s",kqacut?"w/ QA":"w/o QA");

  const TString cutCCpTRec = Form("(fMuonNRec > 0 %s)", kqacut? Form(" && %s",(GetQACut()).Data()): Form(" && %s",(GetPtResCut()).Data()));
  const TString conCCpTRec = Form("#mu^{-} rec %s",kqacut?"w/ p_{t}+PID QA":"w/ p_{t} QA");

  const TString cutCCSim = "( fNeutrinoType==14 && fNeutMode>0 && fNeutMode<30 && fMuonNSim>0 )";
  const TString conCCSim = "#nu_{#mu}CC";

  const TString cutIncl = "1";
  const TString conIncl = "";

  const TString cutCCFull = Form("%s && %s", cutCCSim.Data(), cutCCRec.Data());
  const TString conCCFull = conCCSim+conCCRec;

  const TString cutSimFull = xvar.Contains("Sim")? cutCCSim : cutCCFull;
  const TString conSimFull = xvar.Contains("Sim")? conCCSim : conCCFull;

  TString cuts[10];
  TString cname[10];
  Int_t ncut = 0;

  TString cut0;
  TString header;

  //================================================================================
  if(cvar == "CmpOsci" && pwd.Contains("SK")){
    ncut = 2;

    cuts[0]="fkOsciMuSurvive==1";
    cname[0]="#nu_{#mu}#rightarrow#nu_{#mu}";

    cuts[1]="fkOsciMuSurvive==0";
    cname[1]="#nu_{#mu}#rightarrow others";

    //everything
    cut0 = cutIncl;
    header=conIncl;
  }
  else if(cvar == "EnuNCCC"){
    ncut=2;

    cuts[0]= cutCCSim;
    cname[0]=conCCSim;

    cuts[1]= Form("!(%s)", cuts[0].Data());
    cname[1]=Form("#slash{%s}",cname[0].Data());

    xmin = 5e-2;

    //everything
    cut0 = cutIncl;
    header=conIncl;
  }
  //================================================================================
  else if(cvar == "CCeff"){
    ncut=2;

    cuts[0]= cutCCRec; 
    cname[0]=conCCRec;

    cuts[1]= Form("!(%s)", cuts[0].Data());
    cname[1]=Form("#slash{%s}",cname[0].Data());

    cut0 = cutCCSim;
    header=conCCSim;
  }
  //================================================================================
  else if(cvar == "CCpurity"){
    ncut=2;

    cuts[0]= cutCCSim;
    cname[0]=conCCSim;

    cuts[1]= Form("!(%s)", cuts[0].Data());
    cname[1]=Form("#slash{%s}",cname[0].Data());

    cut0 = cutCCRec;
    header = conCCRec;
  }
  //================================================================================
  else if(cvar == "PID" && pwd.Contains("ND280")){
    Int_t id = 0;

    cuts[id]="fMuonSimPDG==13";
    cname[id++]="#mu^{-}";
    
    cuts[id]="fMuonSimPDG==-13";
    cname[id++]="#mu^{+}";

    cuts[id]="abs(fMuonSimPDG)==11";
    cname[id++]="e";

    cuts[id]="abs(fMuonSimPDG)==211";
    cname[id++]="#pi";
    
    cuts[id]="abs(fMuonSimPDG)==2212";
    cname[id++]="p";
    
    ncut = id;

    //require also CCSim to avoid double counting with CCpurity
    cut0 = cutCCpTRec;
    header = conCCpTRec;
  }
  //================================================================================
  else if(cvar == "PtRes"){
    Int_t id = 0;
  
    cuts[id]="(fMuonRecPt->Mag()/fMuonSimPt->Mag())>=2";
    cname[id++]=">2";
    
    cuts[id]="abs(fMuonRecPt->Mag()/fMuonSimPt->Mag()-1)>0.2 && (fMuonRecPt->Mag()/fMuonSimPt->Mag())<2";
    cname[id++]="<0.8, 1.2-2";
    
    cuts[id]="abs(fMuonRecPt->Mag()/fMuonSimPt->Mag()-1)<=0.2 && abs(fMuonRecPt->Mag()/fMuonSimPt->Mag()-1)>0.1";
    cname[id++]="0.8-0.9, 1.1-1.2";
    
    cuts[id]="abs(fMuonRecPt->Mag()/fMuonSimPt->Mag()-1)<=0.1";
    cname[id++]="0.9-1.1";

    ncut = id;

    cut0= cutCCFull;
    header=conCCFull+" #frac{p_{t}^{#mu rec}}{p_{t}^{#mu sim}}";
  }
  else if(cvar == "PRes"){
    Int_t id = 0;
  
    cuts[id]="(fMuonRec->P()/fMuonSim->P())>=2";
    cname[id++]=">2";
    
    cuts[id]="abs(fMuonRec->P()/fMuonSim->P()-1)>0.2 && (fMuonRec->P()/fMuonSim->P())<2";
    cname[id++]="<0.8, 1.2-2";
    
    cuts[id]="abs(fMuonRec->P()/fMuonSim->P()-1)<=0.2 && abs(fMuonRec->P()/fMuonSim->P()-1)>0.1";
    cname[id++]="0.8-0.9, 1.1-1.2";
    
    cuts[id]="abs(fMuonRec->P()/fMuonSim->P()-1)<=0.1";
    cname[id++]="0.9-1.1";

    ncut = id;

    cut0= cutCCFull;
    header=conCCFull+" #frac{p_{#mu rec}}{p_{#mu sim}}";
  }
  else if(cvar == "MultiplicityFraction"){
    ncut=2;

    cuts[0]="fMultiplicity==1";
    cname[0]="N_{all rec}=1";

    cuts[1]="fMultiplicity>1";
    cname[1]="N_{all rec}>1";

    cut0 = cutCCFull;
    header=conCCFull;
  }
  else if(cvar == "MulPSelQ2" ){
    ncut=3;

    cuts[2]="fMultiplicity==1";
    cname[2]="N_{all rec}=1";

    cuts[1]="fMultiplicity>1 && fMuonRec->P()<=1";
    cname[1]="N_{all rec}>1, p_{#mu rec}#leq 1";

    cuts[0]="fMultiplicity>1 && fMuonRec->P()>1";
    cname[0]="N_{all rec}>1, p_{#mu rec}>1";

    cut0 = cutCCFull;
    header=conCCFull;
  }
  else if(cvar == "PtSelQ2"){
    ncut=2;

    cuts[1]="fMuonRecPt->Mag()<=1";
    cname[1]="p_{t}^{#mu rec}#leq 1";

    cuts[0]="fMuonRecPt->Mag()>1";
    cname[0]="p_{t}^{#mu rec}>1";

    cut0 = cutCCFull;
    header=conCCFull;
  }
  //================================================================================
  else if(cvar == "ModePPt"){
    Int_t id = 0;

    TString modename;

    cuts[id]=GetModeCut(1,modename);
    cname[id++]=modename;

    cuts[id]=GetModeCut(2,modename);
    cname[id++]=modename;

    cuts[id]=GetModeCut(3,modename);
    cname[id++]=modename;

    cuts[id]=GetModeCut(4,modename);
    cname[id++]=modename;

    cuts[id]=GetModeCut(5,modename);
    cname[id++]=modename;

    //checked, no others!
    //cuts[id]=Form("!(%s) && !(%s) && !(%s) && !(%s)", cuts[0].Data(), cuts[1].Data(), cuts[2].Data(), cuts[3].Data());
    //cname[id++]="others";

    ncut=id;
    //for Rec, the NC component has mode >=30
    //for Sim, checked that NC ~ MuonNSim=0

    //cut0 =cutCCRec;
    //header = conCCRec;
    cut0 =cutSimFull;
    header = conSimFull;
  }
  else if(cvar == "ParentPPt" && pwd.Contains("ND280")){
    const Int_t parents[]={13, 130,211,321};
    const TString parname[]={"#mu", "K_{L}","#pi^{#pm}","K^{#pm}"};

    ncut = sizeof(parents)/sizeof(Int_t);

    for(Int_t icut=0; icut<ncut; icut++){
      cuts[icut]=Form("TMath::Abs(fNeutrinoParentPDG)==%d", parents[icut]);
      cname[icut]=parname[icut];
    }

    cut0 =cutSimFull;
    header = conSimFull;
  }
  else if(cvar == "Q2PPt" || cvar == "EnuPPt" ){
    TString tmpcv, tmpcn;
    if(cvar.Contains("Q2")){
      tmpcv = "Q^{2}";
      tmpcn = "fQ2Sim";
    }
    else if(cvar.Contains("Enu")){
      tmpcv = "E_{#nu}";
      tmpcn = "fNeutrinoSim->E()";
    }

    const Double_t nodes[]={-1e10, 0.2, 0.4, 1, 2, 1e10};
    ncut = sizeof(nodes)/sizeof(Double_t)-1;
    for(Int_t icut=0; icut<ncut; icut++){
      TString & cncn = cname[icut];
      TString & ctct = cuts[icut];
      if(fabs(nodes[icut])<10){
        cncn=Form("%.1f<%s", nodes[icut], tmpcv.Data());
      }
      else{
        cncn=Form("%s", tmpcv.Data());
      }

      if(fabs(nodes[icut+1])<10){
        cncn+=Form("<%.1f", nodes[icut+1]);
      }

      ctct=Form("%s>%f && %s<%f", tmpcn.Data(), nodes[icut], tmpcn.Data(), nodes[icut+1]);
    }


    cut0 =cutSimFull;
    header = conSimFull;
  }
  else{
    printf("wrong cvar %s\n", cvar.Data()); exit(1);
  }

  //add survival requirement, i.e. physics case at SK
  if(kosci && cvar!="CmpOsci"){
    //there is osci = -999 for ND
    cut0 += "&& fkOsciMuSurvive==1 ";
  }

  //===>
  if(xvar=="fNeutrinoSim->E()"){
    xmin = 0.99e-1;
  }
  
  if(xvar=="fMuonRec->P()"){
    //xmin = 0.99e-1;
    xmin = 1e-3;
  }

  if(xvar=="fMultiplicity"){
    xmin = 0.5;
    
    if(pwd.Contains("SK")){
      xmax = 5.5;
    }
    else{
      xmax=10.5;
    }
    xnbin = xmax-xmin;
  }
  
  if(xvar=="fNeutMode"){
    xmin = -61.5;
    xmax =  61.5;
    xnbin = xmax-xmin;
  }

  if(xvar=="fNeutrinoType"){
    xmin = -15.5;
    xmax =  15.5;
    xnbin = xmax-xmin;
  }

  if(xvar=="(fMuonCharge[1]==-1)*2+(fMuonCharge[0]==-1)"){
    xmin = -0.5;
    xmax =  3.5;
    xnbin = xmax-xmin;
  }

  if(xvar=="fMuonRecFlightPath->Pt()"){
    xmin = 0;
    xmax = 5;
    xnbin = 25;
  }

  if(xvar=="fMuonRecVertex->X()" || xvar=="fMuonRecVertex->Y()"){
    xmin = -1;
    xmax = 1;
    xnbin = 50;
  }

  if(xvar=="fMuonRecVertex->Z()"){
    xmin = 0.1;
    xmax = 0.5;
    xnbin = 50;
  }

  if(xvar=="fMuonRecEndPos->X()" || xvar=="fMuonRecEndPos->Y()"){
    xmin = -3;
    xmax = 3;
    xnbin = 60;
  }

  if(xvar=="fMuonRecEndPos->Z()"){
    xmin = -0.5;
    xmax = 4.5;
    xnbin = 80;
  }

  if(xvar=="fMuonRecFlightPath->Mag()"){
    xmin = 0;
    xmax = 6;
    xnbin = 30;
  }


  if(xvar=="fRangeMomentumMuon"){
    xmin = 0;
    xmax = 2;
    xnbin = 30;
  }

  if(xvar=="fMuonRecNhits[0]"){
    xmin = 0;
    xmax = 300;
    xnbin=150;
  }

  if(xvar=="fMuonLength[0]"||xvar=="fMuonLength[1]"||xvar=="fMuonLength[2]"){
    xmin = -0.1;
    xmax = 3.1;
    xnbin=150;
  }

  if(xvar=="fMuonRecNhits[1]"){
    xmin = 0;
    xmax = 150;
    xnbin=75;
  }

  if(xvar=="fMuonTypeSim"){
    xmin = -6.5;
    xmax = 6.5;
    xnbin = xmax-xmin;
  }

  if(xvar=="fNeutrinoParentType"){
    xmin = -6.5;
    xmax = 6.5;
    xnbin = xmax-xmin;
  }

  if(xvar=="fNeutrinoRec->Theta()*TMath::RadToDeg()"){
    xmin = 1.3;
    xmax = 2.0;
    xnbin = 50;
  }

  if(xvar=="fNeutrinoSim->Theta()*TMath::RadToDeg()"){
    xmin = 0;
    xmax = 3;
    xnbin = 50;
  }

  if(xvar=="fMuonAlphaSim" || xvar=="fMuonAlphaRec" || xvar=="fMuonRec->Theta()*TMath::RadToDeg()"){
    xmin = 0;
    xmax = 180;
    xnbin = 50;
  }

  if((cvar=="ModePPt" || cvar=="Q2PPt") ){
    if(xvar=="fMuonSimPt->Mag()"){
      xmin = 8e-2;
      xmax = 3;
      xnbin = 30;
    }
    else if(xvar=="fMuonSim->P()"){
      xmin = 8e-2;
      xmax = 10;
      xnbin = 15;
    }
  }

  //<===

  //check!
  //should only be a very small amount (~10) underflow
  //grep FLOW tmp | grep -v "UNDERFLOW 0 "  | grep UNDERFLOW
  //the non-0 is due to MuonN cut, and should all be the same value for the same cut: 3 possibilities: fMuonNSim>0, fMuonNSim>0 && fMuonNRec>0, fMuonNRec>0.
  //grep last tmp  | grep "d 0;" -v| grep " d "

  //----------------------------------------------------------------------
  const Int_t cols[]={kRed-4, kOrange, kCyan, kMagenta-7, kGreen};
  const Int_t mks[]={20,21,22,23,29};
  const TString tottit(Form(";%s;Counts", (getTitleFromVar(xvar)).Data()));
  TLegend *lg = new TLegend(0.67, 0.72, 0.95, 0.90);
  style::ResetStyle(lg);
  lg->SetHeader(header);
  lg->SetBorderSize(0);
  lg->SetFillColor(0);

  TH1D *hcount[]={0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0};
  //only for calculation of efficiency
  TH1D * htot = new TH1D("htot","", xnbin, xmin, xmax); 
  if(xvar!="fMultiplicity"){
    NeutrinoTools::BinLog(htot->GetXaxis());
  }
  style::ResetStyle(htot);

  THStack * countStack = new THStack("countStack","");
  countStack->SetTitle(tottit);
  
  Double_t totn=0;
  for(Int_t ii=0; ii<ncut; ii++){
    printf("[%d]: cuts %s , cname %s\n", ii, cuts[ii].Data(), cname[ii].Data());

    //some fMuonSimPt->Mag has ~ 1e-10 entries (~4), all with fNeutMode = 26, very parallel to neutrino
    
    hcount[ii] = new TH1D(Form("hh%d",ii),"", xnbin, xmin, xmax); 
    if(xvar!="fMultiplicity"){
      NeutrinoTools::BinLog(hcount[ii]->GetXaxis());
    }
    style::ResetStyle(hcount[ii]);

    const TString finalcut=(cut0+" && ")+cuts[ii];
    tree->Draw(xvar+Form(">>hh%d",ii), finalcut);
   
    totn+=hcount[ii]->GetEntries();
    printf("\nxvar %s || cvar %s || finalcut { %s } || cuts[%d] %s || kqacut %d || %s %d: entries %.0f/%.0f d %.0f; UNDERFLOW %.0f OVERFLOW %.0f\n\n",xvar.Data(), cvar.Data(), finalcut.Data(), ii, cuts[ii].Data(), kqacut,  ii==ncut-1? "last":"", ii, hcount[ii]->GetEntries(), tree->GetEntries(), totn-tree->GetEntries(), hcount[ii]->GetBinContent(0), hcount[ii]->GetBinContent(hcount[ii]->GetNbinsX()+1));

    hcount[ii]->SetFillColor(cols[(ii+1)%5]);
    hcount[ii]->SetLineColor(kBlack);
    //hcount[ii]->SetLineColor(cols[(ii+1)%5]);
    hcount[ii]->SetMarkerColor(cols[(ii+1)%5]); 
    hcount[ii]->SetMarkerStyle(mks[ii]);
    hcount[ii]->SetMarkerSize(2);
    hcount[ii]->SetLineWidth(1);

    hcount[ii]->SetTitle(cname[ii]+tottit);
    countStack->Add(hcount[ii]);
    lg->AddEntry(hcount[ii], cname[ii],"f");

    htot->Add(hcount[ii]);
  }

  //tmp
  countStack->Draw();
  //has to draw first
  style::ResetStyle(countStack);

  TString effTit(tottit);
  //effTit.ReplaceAll("Counts","#epsilon");
  effTit.ReplaceAll("Counts","Fraction");
  THStack * effStack = new THStack("effStack","");
  effStack->SetTitle(effTit);
  htot->Sumw2();
  for(Int_t ii=0; ii<ncut; ii++){
    printf("doing efficiency [%d]: %s , %s\n", ii, cuts[ii].Data(), cname[ii].Data());
    TH1D * heff = (TH1D*) hcount[ii]->Clone(Form("%seff", hcount[ii]->GetName()));
    style::ResetStyle(heff);
    heff->Sumw2();
    heff->Divide(heff, htot, 1, 1, "B");
    heff->SetTitle(cname[ii]+effTit);
    effStack->Add(heff);
  }

  //-------------------------------------------------------------
  
  TString tmp(xvar);
  tmp.ReplaceAll("*TMath::RadToDeg()","");
  tmp.ReplaceAll(">","");
  tmp.ReplaceAll("(","");
  tmp.ReplaceAll(")","");
  tmp.ReplaceAll("-","");
  tmp.ReplaceAll("/","");
  tmp.ReplaceAll("[","");
  tmp.ReplaceAll("]","");
  tmp.ReplaceAll("+","");
  tmp.ReplaceAll("*","");
  tmp.ReplaceAll("=","");
  tmp.Prepend("outplot/");
  printf("tmp %s\n", tmp.Data());

  //c1->SetTopMargin(0.08);
 
  for(int ii=0; ii<2; ii++){
    gPad->SetLogy(ii); 
    countStack->SetMinimum(ii?0.9:0);
    for(Int_t jj=0; jj<2; jj++){
      gPad->SetLogx(jj);
      countStack->Draw("hist");
      lg->Draw();
      c1->Print(tmp+Form("countWithStack%sLogx%dLogy%dqacut%d%s.png",cvar.Data(), gPad->GetLogx(),  gPad->GetLogy(), kqacut,kosci?"MuSurv":""));
      
      countStack->Draw("nostack lpe");
      lg->Draw();
      c1->Print(tmp+Form("countNoStack%sLogx%dLogy%dqacut%d%s.png", cvar.Data(), gPad->GetLogx(),  gPad->GetLogy(), kqacut,kosci?"MuSurv":""));
    }
  }

  //tmp
  effStack->Draw();
  style::ResetStyle(effStack);

  //c1->SetTopMargin(0.02);

  for(int ii=0; ii<2; ii++){
    gPad->SetLogy(ii); 
    //effective, but the number is not accurate
    effStack->SetMinimum(ii?1e-3:0);
    effStack->SetMaximum(ii?15:1.41);
  
    for(Int_t jj=0; jj<2; jj++){
      gPad->SetLogx(jj);
      effStack->Draw("hist");
      lg->Draw();
      c1->Print(tmp+Form("effWithStack%sLogx%dLogy%dqacut%d%s.png",cvar.Data(), gPad->GetLogx(),  gPad->GetLogy(), kqacut,kosci?"MuSurv":""));
      
      effStack->Draw("nostack lpe");
      lg->Draw();
      c1->Print(tmp+Form("effNoStack%sLogx%dLogy%dqacut%d%s.png", cvar.Data(), gPad->GetLogx(),  gPad->GetLogy(), kqacut,kosci?"MuSurv":""));
    }
  }
}

void plotPtInvRes(const TString formula, const Int_t kInv=kTRUE, 
                  const Int_t nhitcut=0, const Double_t larm = 0, const Double_t maxp = 0, const Double_t maxchi = 0,
                  const Int_t kCharge=0, const Double_t fpath = 0, const Int_t maxnhit = 0, const Double_t minrange = 0,
                  const Bool_t kZaxis=kFALSE, const Bool_t kTrueMu = kFALSE, const Int_t kNuDirOpt = 0
                  )
{
  const TString pwd=gSystem->pwd();

  if(pwd.Contains("SK")){
    if(kZaxis || kTrueMu || kNuDirOpt!=1 ){
      return;
    }

    if(nhitcut || larm>0 || maxp>0 || maxchi>0 || kCharge>0 || fpath>0 || maxnhit>0 || minrange>0){
      return;
    }      
  }
  else{
    if(kZaxis){
      //require PID, others not care
      if(kTrueMu==0){
        return;
      }
    }

    if(kTrueMu){
      if(kNuDirOpt!=2){
        return;
      }
    }
  }

  const TString tag0=Form("kZaxis%d_kTrueMu%d_kNuDirOpt%d", kZaxis, kTrueMu, kNuDirOpt);
  const TString tag=Form("nhit%02d_larm%03.0f_maxp%03.0f_maxchi%03.0f_kCharge%d_fpath%03.0f_maxnhit%d_minrange%03.0f_formula%s_kInv%d_%s",nhitcut, larm*100, maxp*100, maxchi*100, kCharge, fpath*100, maxnhit, minrange*100, formula.Data(),  kInv, tag0.Data());

  //================================================

  gStyle->SetOptStat(0);
  
  gROOT->LoadMacro("style.cxx+");

  style::SetGlobalStyle();

  style::fgkYTitleOffset = 1.4;
  style::fgkXTitleOffset = 1.4;

  TCanvas *c0= new TCanvas;
  style::PadSetup(c0);

  const Double_t lm=0.15;
  Double_t rm=0.13;
  const Double_t tm=0.13;
  const Double_t bm=0.17;

  c0->SetLeftMargin(lm);
  c0->SetRightMargin(rm);
  c0->SetTopMargin(tm);
  c0->SetBottomMargin(bm);

  gStyle->SetOptTitle(1);
  //gStyle->SetTitleX(0.05);
  gStyle->SetTitleW(0.99);

  const Int_t xnbin = 20;
  const Double_t xmin = 0.05;
  const Double_t xmax = kNuDirOpt==3? 3.5: 2;

  const Int_t ynbin = 120*(kInv?1:2);
  //const Double_t ymax = (kInv? 0.6 : 0.3)*(pwd.Contains("SK")?1:2);
  const Double_t ymax = 1.2;//(kInv? 0.6 : 0.3)*2;
  const Double_t ymin = -ymax;
  //TH2D * h2=new TH2D("h2",cut, pwd.Contains("SK")?50:10,  pwd.Contains("SK")?0.1:0.2, 2, ynbin, ymin, ymax); NeutrinoTools::BinLog(h2->GetXaxis());
  TH2D * h2=new TH2D("h2","", xnbin, xmin, xmax, ynbin, ymin, ymax); NeutrinoTools::BinLog(h2->GetXaxis());

  TH1D * halleff = new TH1D("halleff","",xnbin, xmin, xmax); NeutrinoTools::BinLog(halleff->GetXaxis());
  TH1D * heff = new TH1D("heff","",xnbin, xmin, xmax); NeutrinoTools::BinLog(heff->GetXaxis());

  TH1D * hallmuonpurity = new TH1D("hallmuonpurity","",xnbin, xmin, xmax); NeutrinoTools::BinLog(hallmuonpurity->GetXaxis());
  TH1D * hmuonpurity = new TH1D("hmuonpurity","",xnbin, xmin, xmax); NeutrinoTools::BinLog(hmuonpurity->GetXaxis());

  TH1D * hallnupurity = new TH1D("hallnupurity","",xnbin, xmin, xmax); NeutrinoTools::BinLog(hallnupurity->GetXaxis());
  TH1D * hnupurity = new TH1D("hnupurity","",xnbin, xmin, xmax); NeutrinoTools::BinLog(hnupurity->GetXaxis());

  //================================================

  const TString cut0("fMuonNRec>0");
  const TString cut= cut0
    + (nhitcut>0? Form(" && fMuonRecNhits[0]>%d", nhitcut):"") 
    + (larm>0?    Form(" && fMuonRecFlightPath->Pt()>%.2f", larm):"")
    + (maxp>0?    Form(" && fMuonRec->P()<%.2f",maxp):"")
    + (maxchi>0?  Form(" && %.2f<(fMuonRecChi2/fMuonRecNDOF) && (fMuonRecChi2/fMuonRecNDOF)<%.2f", 1/maxchi, maxchi):"")
    + (kCharge>0? Form(" && ((fMuonCharge[1]==-1)*2+(fMuonCharge[0]==-1))==%d", kCharge):"")
    + (fpath>0?   Form(" && fMuonRecFlightPath->Mag()>%.2f", fpath):"")
    + (maxnhit>0? Form(" && fMuonRecNhits[0]<%d", maxnhit):"")
    + (minrange>0?Form(" && fRangeMomentumMuon>%.2f", minrange):"")
    + (kTrueMu ? " && fMuonSimPDG==13" : "")
    + ( (kNuDirOpt==3) ? "&& (fNeutMode==1||fNeutMode==2)" : "" ); //momentum resolution only relevant for CCQE, for Enu resolution


  TString xRec, xSim;

  if(kNuDirOpt==0){  
    //pt w.r.t. (001) 
    xRec = "fMuonRec->Pt()";
  }
  else if(kNuDirOpt==1){
    //using parent dec point - vertex
    xRec = "fMuonRecPt->Mag()";
  }
  else if(kNuDirOpt==2){
    //using true nutrino direction
    xRec = "fMuonPtGuess->Mag()";
  }
  else if(kNuDirOpt!=3 && kNuDirOpt!=4){
    printf("\n!!wrong kNuDirOpt %d\n\n", kNuDirOpt);exit(1);
  }
  //true pt
  xSim = "fMuonSimPt->Mag()";

  if(kNuDirOpt==3){
    xRec = "fMuonRec->P()";
    xSim = "fMuonSim->P()";
  }

  //overwrite kNuDirOpt
  if(kZaxis){
    xRec = "fMuonRec->Pt()";
    xSim = "fMuonSim->Pt()";
  }

  //overwrite kNuDirOpt again
  if(kNuDirOpt==4){
    //testtest the displayed x title is still not correct!
    xRec = "fDeltaRecPt->Mag()";
    xSim = "fDeltaSimPt->Mag()";
  }

  TString plvar=Form("%s:%s>>h2",kInv? "1/rec-1/sim":"rec/sim-1","sim");
  plvar.ReplaceAll("rec",xRec);
  plvar.ReplaceAll("sim",xSim);

  if(kNuDirOpt==4){
    plvar.ReplaceAll(":fDelta",":fMuon");
  }

  printf("check tag %s || plvar %s || cut { %s }\n\n",tag. Data(), plvar.Data(), cut.Data());
  tree->Draw(plvar,cut+"&& fMuonNSim>0");

  tree->Draw(xSim+">>halleff");
  tree->Draw(xSim+">>heff",cut);

  tree->Draw(xRec+">>hallmuonpurity", cut);
  tree->Draw(xRec+">>hmuonpurity", cut+" && fMuonSimPDG==13");

  tree->Draw(xRec+">>hallnupurity", cut);
  tree->Draw(xRec+">>hnupurity", cut+" && fNeutrinoType==14 && fNeutMode<30 && fNeutMode>0");

  const TString xTitRec = getTitleFromVar(xRec);
  const TString xTitSim = getTitleFromVar(xSim);
  //============================================================================
  
  TString atit(cut);
  atit.ReplaceAll("((fMuonCharge[1]==-1)*2+(fMuonCharge[0]==-1))==1","tpccharge01(-+)");
  atit.ReplaceAll("((fMuonCharge[1]==-1)*2+(fMuonCharge[0]==-1))==2","tpccharge01(+-)");
  atit.ReplaceAll("((fMuonCharge[1]==-1)*2+(fMuonCharge[0]==-1))==3","tpccharge01(--)");
  atit.ReplaceAll("&&"," ");
  atit.ReplaceAll("   "," ");
  atit.ReplaceAll("  "," ");
  atit.ReplaceAll("fMuonNSim>0","");
  atit.ReplaceAll("fMuonNRec>0","");
  atit.ReplaceAll("fMuonRecNhits[0]","tpcnhit0");
  atit.ReplaceAll("fMuonRecNhits[1]","tpcnhit1");
  atit.ReplaceAll("fMuonRecFlightPath->Pt()","larm");
  atit.ReplaceAll("fMuonCharge[0]","tpccharge0");
  atit.ReplaceAll("fMuonCharge[1]","tpccharge1");
  atit.ReplaceAll("fMuonSimPDG==13","#mu^{-}");
  atit.ReplaceAll("fMuonRec->P()","p_{#mu}");
  atit.ReplaceAll("fMuonRecFlightPath->Mag()","fpath");
  atit.ReplaceAll("(fNeutMode==1||fNeutMode==2)","CCQE");
  atit.ReplaceAll("(fMuonRecChi2/fMuonRecNDOF) (fMuonRecChi2/fMuonRecNDOF)","#chi^{2}_{track}/NDOF");
  atit.ReplaceAll("fRangeMomentumMuon","R_{#mu rec}");

  TString tottitle(atit+Form(" %s;%s;%s", tag0.Data(), xTitSim.Data(), kInv?"1/rrr-1/sss (c/GeV)":"rrr/sss-1"));
  tottitle.ReplaceAll("sss", "p_{t}^{#mu sim}");
  tottitle.ReplaceAll("rrr", "p_{t}^{#mu rec}");

  tottitle.ReplaceAll("kZaxis0","");
  tottitle.ReplaceAll("kZaxis1_kTrueMu1_kNuDirOpt2","intrinsic p_{t}");

  tottitle.ReplaceAll("_k","k");

  tottitle.ReplaceAll("kTrueMu0","");
  if(tottitle.Contains("#mu^{-}")){
    tottitle.ReplaceAll("kTrueMu1","");
  }

  tottitle.ReplaceAll("kNuDirOpt0","#vec{#nu}=(001)");
  tottitle.ReplaceAll("kNuDirOpt1","#vec{#nu}=D-V");
  tottitle.ReplaceAll("kNuDirOpt2","#vec{#nu}_{sim}");
  tottitle.ReplaceAll("kNuDirOpt3","");


  h2->SetTitle(tottitle);

  //printf("check overflow and underflow: %f %f\n", h2->Integral(0,h2->GetNbinsX(), 0, 0), h2->Integral(0,h2->GetNbinsX(), h2->GetNbinsY(), h2->GetNbinsY()));

  h2nor=NeutrinoTools::NormalHist(h2, 5, 1);
  gPad->SetLogx();
  style::ResetStyle(h2nor);
  h2nor->Draw("colz");
  
  c0->Print(Form("outplot/ptinvres2d%s.png", tag.Data()));

  gPad->SetLogz();
  c0->Print(Form("outplot/ptinvres2d%slogz.png", tag.Data()));

  //===========================
  const Double_t thres = 50;
  TList * lout=new TList;
  lout->Add(h2);

  TH1D *hrmsnor=0x0, *hrmsmpv=0x0, *hrmswid=0x0, *hrmsres=0x0, *hrmschi=0x0;
  NeutrinoTools::FitSlicesY(h2, hrmsnor, hrmsmpv, hrmswid, hrmsres, hrmschi, "RMS", thres, lout);

  TH1D *hnor=0x0, *hmpv=0x0, *hwid=0x0, *hres=0x0, *hchi=0x0;
  NeutrinoTools::FitSlicesY(h2, hnor, hmpv, hwid, hres, hchi, formula, thres, lout);

//  h2->FitSlicesY();
//  hmpv=h2_1;
//  hwid=h2_2;
  
  TFile * fout=new TFile(Form("outplot/fitslice%s.root", tag.Data()),"recreate");
  lout->Write();
  fout->Save();
  fout->Close();
  delete fout;
  //============================================================================
  rm = 0.02;

  TCanvas *c1=new TCanvas;
  style::PadSetup(c1);

  gPad->SetLogx();

  c1->SetLeftMargin(lm);
  c1->SetRightMargin(rm);
  c1->SetTopMargin(tm);
  c1->SetBottomMargin(bm);

  gPad->SetGrid();
  hmpv->SetTitle(tottitle);
  hmpv->SetYTitle(Form("<%s>", h2->GetYaxis()->GetTitle()));
  hmpv->SetMinimum(-0.2);
  hmpv->SetMaximum(0.2);
  style::ResetStyle(hmpv);

  hmpv->SetMarkerStyle(20);
  hmpv->SetMarkerSize(1);
  hmpv->SetMarkerColor(kBlack);
  hmpv->SetLineWidth(1);
  hmpv->SetLineColor(kBlack);
  hmpv->Draw("e");
  c1->Print(Form("outplot/ptinvresmean%s.png", tag.Data()));

  //============================================================================

  TCanvas *c2 =new TCanvas;
  style::PadSetup(c2);

  c2->SetLeftMargin(lm);
  c2->SetRightMargin(rm);
  c2->SetTopMargin(tm);
  c2->SetBottomMargin(bm);

  gPad->SetLogx(1);
  gPad->SetGrid();

  hwid->SetTitle(tottitle);
  hwid->SetYTitle(Form("#sigma(%s)", h2->GetYaxis()->GetTitle()));
  hwid->SetMinimum(0);
  //hwid->SetMaximum(pwd.Contains("SK")?0.2:0.4);//0.6
  hwid->SetMaximum(kInv?0.4:0.2);

  style::ResetStyle(hwid);
  style::ResetStyle(hrmswid);
  
  hwid->SetMarkerStyle(20);
  hwid->SetMarkerSize(1);
  hwid->SetMarkerColor(kBlack);
  hwid->SetLineWidth(1);
  hwid->SetLineColor(kBlack);

  hrmswid->SetMarkerStyle(22);
  hrmswid->SetMarkerSize(1.5);
  hrmswid->SetMarkerColor(kRed);
  hrmswid->SetLineColor(kRed);
  hrmswid->SetLineWidth(1);

  hwid->Draw("e");
  hrmswid->Draw("same e");

  TLegend * wl = new TLegend(0.15,0.7,0.4,0.87);
  style::ResetStyle(wl,0.2);
  wl->SetTextAlign(12);
  //wl->SetBorderSize(1);
  wl->AddEntry(hrmswid, "RMS","lp");
  wl->AddEntry(hwid,formula,"lp");
  wl->Draw();

  c2->Print(Form("outplot/ptinvreswid%s.png", tag.Data()));

  if(kInv!=1 || (kNuDirOpt!=1 && kNuDirOpt!=3 && kNuDirOpt!=4) || kZaxis!=0 ||  kTrueMu!=0){
    return;
  }
  //============================================================================

  TCanvas *c3 =new TCanvas;
  style::PadSetup(c3);

  c3->SetLeftMargin(lm);
  c3->SetRightMargin(rm);
  c3->SetTopMargin(tm);
  c3->SetBottomMargin(bm);

  gPad->SetLogx(1);
  gPad->SetGrid();

  halleff->Sumw2();
  heff->Sumw2();
  heff->Divide(heff, halleff,1,1,"B");

  heff->SetTitle(Form(" %s;%s; efficiency",hwid->GetTitle(), xTitSim.Data()));
  heff->SetMinimum(gPad->GetLogy()?0.01:0);
  heff->SetMaximum(heff->GetBinContent(heff->GetMaximumBin())<0.4?0.5:1.1);
  style::ResetStyle(heff);

  heff->SetMarkerStyle(20);
  heff->SetMarkerSize(1);
  heff->SetMarkerColor(kBlack);
  heff->SetLineWidth(1);
  heff->SetLineColor(kBlack);
  heff->Draw("e");
  c3->Print(Form("outplot/ptinvreseff%s.png", tag.Data()));

  //============================================================================

  TCanvas *c4 =new TCanvas;
  style::PadSetup(c4);

  c4->SetLeftMargin(lm);
  c4->SetRightMargin(rm);
  c4->SetTopMargin(tm);
  c4->SetBottomMargin(bm);

  gPad->SetLogx(1);
  gPad->SetGrid();

  hallmuonpurity->Sumw2();
  hmuonpurity->Sumw2();
  hmuonpurity->Divide(hmuonpurity, hallmuonpurity,1,1,"B");

  //it is rec!!
  hmuonpurity->SetTitle(Form(" %s;%s; #mu^{-} purity",hwid->GetTitle(), xTitRec.Data()));
  hmuonpurity->SetMinimum(0.7);
  hmuonpurity->SetMaximum(1.1);
  style::ResetStyle(hmuonpurity);

  hmuonpurity->SetMarkerStyle(20);
  hmuonpurity->SetMarkerSize(1);
  hmuonpurity->SetMarkerColor(kBlack);
  hmuonpurity->SetLineWidth(1);
  hmuonpurity->SetLineColor(kBlack);
  hmuonpurity->Draw("e");
  c4->Print(Form("outplot/ptinvresmuonpurity%s.png", tag.Data()));

  //============================================================================

  TCanvas *c5 =new TCanvas;
  style::PadSetup(c5);

  c5->SetLeftMargin(lm);
  c5->SetRightMargin(rm);
  c5->SetTopMargin(tm);
  c5->SetBottomMargin(bm);

  gPad->SetLogx(1);
  gPad->SetGrid();

  hallnupurity->Sumw2();
  hnupurity->Sumw2();
  hnupurity->Divide(hnupurity, hallnupurity,1,1,"B");

  //it is rec!!
  hnupurity->SetTitle(Form(" %s;%s; #nu_{#mu} CC purity",hwid->GetTitle(), xTitRec.Data()));
  hnupurity->SetMinimum(0.7);
  hnupurity->SetMaximum(1.1);
  style::ResetStyle(hnupurity);

  hnupurity->SetMarkerStyle(20);
  hnupurity->SetMarkerSize(1);
  hnupurity->SetMarkerColor(kBlack);
  hnupurity->SetLineWidth(1);
  hnupurity->SetLineColor(kBlack);
  hnupurity->Draw("e");
  c5->Print(Form("outplot/ptinvresnupurity%s.png", tag.Data()));
}

void ptSpectrum(const Bool_t kqacut, const Int_t imode=0)
{
  const TString pwd=gSystem->pwd();

  if(pwd.Contains("SK")){
    if(kqacut){
      return;
    }
  }

  //===============================================

  gStyle->SetOptStat(0);
  
  gROOT->LoadMacro("style.cxx+");

  style::SetGlobalStyle();

  style::fgkYTitleOffset = 1.4;
  style::fgkXTitleOffset = 1.4;

  TCanvas *c0= new TCanvas;
  style::PadSetup(c0);

  gPad->SetLogy();
  gPad->SetLogx();

  Double_t lm=0.15;
  Double_t rm=0.02;
  Double_t tm=0.04;
  Double_t bm=0.17;

  c0->SetLeftMargin(lm);
  c0->SetRightMargin(rm);
  c0->SetTopMargin(tm);
  c0->SetBottomMargin(bm);

  gStyle->SetOptTitle(1);

  const Int_t nxbins = 15;
  const Double_t xmin = 0.05;
  const Double_t xmax = 2;

  //===============================================

  //========>
  const TString xRecVar = kqacut?"fMuonRecPt->Mag()":"fMuonRec->Pt()";

  TString modename;
  const TString modecut = GetModeCut(imode, modename);
  const TString finalcut = Form("(fMuonNRec > 0 %s && %s)", kqacut? Form(" && %s",(GetQACut()).Data()): "", modecut.Data());
  printf("xRecVar %s || finalcut { %s } || kqacut %d \n",xRecVar.Data(), finalcut.Data(), kqacut);
  //<=========

  TH1D * hrec = new TH1D("hrec","",nxbins, xmin, xmax); style::ResetStyle(hrec);   hrec->SetTitle(Form("%s;p_{t}^{#mu} (GeV/c);Counts",""));   
  NeutrinoTools::BinLog(hrec->GetXaxis()); 
  tree->Draw(xRecVar+">>hrec",finalcut);

  TH1D * hsim = new TH1D("hsim","",nxbins, xmin, xmax); style::ResetStyle(hsim); 
  NeutrinoTools::BinLog(hsim->GetXaxis()); 
  tree->Draw("fMuonSimPt->Mag()>>hsim",finalcut);

  TLegend * lg = new TLegend(0.17,0.8,0.25,0.95);
  style::ResetStyle(lg,0.3);
  lg->SetTextAlign(12);
  //lg->SetBorderSize(1);

  lg->SetHeader(Form("%s %s", kqacut?"D-V w/ QA":"#vec{#nu}=(001) w/o QA", modename.Data()));
  hrec->SetMinimum(TMath::Max(1.0,hsim->GetMinimum()*0.5));
  hrec->SetMaximum(hrec->GetBinContent(hrec->GetMaximumBin())*2);
  hrec->SetLineWidth(2);
  hrec->SetLineColor(kBlack);
  hrec->Draw("hist e"); 
  lg->AddEntry(hrec,"rec","l");

  hsim->SetLineWidth(3);
  hsim->SetLineColor(kRed);
  hsim->SetLineStyle(kDashed);
  hsim->Draw("hist e  same");
  lg->AddEntry(hsim,"sim","l");

  lg->Draw();
  TString tag=Form("kqacut%d_mode%s",kqacut, modename.Data());
  tag.ReplaceAll("#","");
  tag.ReplaceAll("{","");
  tag.ReplaceAll("}","");
  tag.ReplaceAll(" or ","");

  c0->Print(Form("outplot/ptSpec_%s.png",tag.Data()));


  c1->SetTopMargin(tm);
  gPad->SetLogx();
  gPad->SetLogy(0);

  gPad->SetGrid();
 
  hrec->Sumw2();
  hsim->Sumw2();
  hrec->Divide(hsim);
  hrec->SetYTitle("spectrum ratio (rec/sim)");
  hrec->SetMarkerStyle(20);
  hrec->SetMarkerSize(1);
  hrec->SetMinimum(0.5);
  hrec->SetMaximum(1.5);
  hrec->Draw("e");

  TLatex *lt=new TLatex(0.55,0.9, lg->GetHeader());
  style::ResetStyle(lt);
  //lt->SetTextAlign(33);
  lt->Draw();

  TF1 * f1=new TF1("f1","1",0,10);
  style::ResetStyle(f1);
  f1->SetLineColor(kRed);
  f1->SetLineWidth(2);
  f1->Draw("same");
  c0->Print(Form("outplot/ptSpRatio_%s.png",tag.Data()));

}

TF1* tmpplotFuncblock(TH1D *hh, TF1 *ftmp, const Int_t kcol, const Int_t lsty,  const TString fitopt, TLegend *lg)
{
  printf("\nfitting %s %s...\n", ftmp->GetName(), ftmp->GetTitle());

  style::ResetStyle(ftmp); 
  ftmp->SetParameters(100,0,0.1); 
  ftmp->SetLineColor(kcol); 
  ftmp->SetLineWidth(3);
  ftmp->SetLineStyle(lsty);
  hh->Fit(ftmp,fitopt); 
  lg->AddEntry(ftmp,ftmp->GetName(),"l");

  return ftmp;
}

plotFunc(const TString fin)
{
  new TFile(fin);

  TH2D *hori = h2;

  const Int_t ibin = 13;
  TH1D *hh = h2->ProjectionY("proj",ibin, ibin);
  const TString header(Form("%s %.2f - %.2f", h2->GetXaxis()->GetTitle(), h2->GetXaxis()->GetBinLowEdge(ibin), h2->GetXaxis()->GetBinUpEdge(ibin)));
  hh->SetTitle(h2_RMSmpv->GetTitle());
  hh->SetYTitle("Counts");
  hh->SetXTitle(hori->GetYaxis()->GetTitle());

  //===

  gROOT->LoadMacro("style.cxx+");

  style::SetGlobalStyle();

  style::fgkYTitleOffset = 1.4;
  style::fgkXTitleOffset = 1.4;

  TCanvas *c0= new TCanvas;
  style::PadSetup(c0);

  Double_t lm=0.15;
  Double_t rm=0.02;
  Double_t tm=0.13;
  Double_t bm=0.17;

  c0->SetLeftMargin(lm);
  c0->SetRightMargin(rm);
  c0->SetTopMargin(tm);
  c0->SetBottomMargin(bm);

  //===

  gStyle->SetOptTitle(1);
  gStyle->SetTitleW(0.99);
  style::ResetStyle(hh);
  hh->SetLineColor(kBlack);
  hh->SetMarkerColor(kBlack);
  hh->SetMarkerStyle(20);
  hh->SetMarkerSize(0.7);
  //gPad->SetLogy();
  hh->SetMaximum(hh->GetBinContent(hh->GetMaximumBin())*(gPad->GetLogy()?3:1.3));

  TLegend * lg = new TLegend(0.17,0.44,0.85,0.83);
  style::ResetStyle(lg,0.07);
  lg->SetTextAlign(12);
  //lg->SetBorderSize(1);
  
  //fits are affected by SetRangeUser!!
  TF1 * fgausML =   tmpplotFuncblock(hh, new TF1("GausML","[0]*TMath::Gaus(x,[1],[2],1)",-2,2),       kRed,     kDashed,     "NL", lg);
  TF1 * fgausLS =   tmpplotFuncblock(hh, new TF1("GausLS","[0]*TMath::Gaus(x,[1],[2],1)",-2,2),       kBlue,    kDotted,     "N",  lg);
  TF1 * fcauchyML = tmpplotFuncblock(hh, new TF1("CauchyML","[0]*TMath::CauchyDist(x,[1],[2])",-2,2), kGreen,   kSolid,      "NL", lg);
  TF1 * fcauchyLS = tmpplotFuncblock(hh, new TF1("CauchyLS","[0]*TMath::CauchyDist(x,[1],[2])",-2,2), kMagenta, kDashDotted, "N",  lg);

  //GetRMS is affected by SetRangeUser
  lg->SetHeader(Form("#splitline{RMS=%.2f, #sigma_{G ML}=%.2f, #sigma_{G LS}=%.2f, #sigma_{C ML}=%.2f, #sigma_{C LS}=%.2f}{%s}", hh->GetRMS(),  fgausML->GetParameter(2), fgausLS->GetParameter(2), fcauchyML->GetParameter(2), fcauchyLS->GetParameter(2),header.Data()));

  //======================================================

  //SetRangeUser only after getting all numbers
  hh->GetXaxis()->SetRangeUser(-0.65,0.65);

  hh->Draw("e");

  TLine * ln = new TLine(0,0,0,hh->GetMaximum());
  ln->SetLineColor(kGray);
  ln->SetLineWidth(1);
  ln->Draw("same");

  fgausML->Draw("same");
  fgausLS->Draw("same");
  fcauchyML->Draw("same");
  fcauchyLS->Draw("same");
  hh->Draw("same");



  lg->Draw();

  TString outn(fin);
  outn=outn(outn.First("/")+1,outn.Length());
  outn=outn(0,outn.First("."));
  cout<<outn<<endl;
  c0->Print(Form("outplot/momDescription%s%s.png", outn.Data(), hh->GetName()));

  //======================================================

  gPad->SetLogy(0);
  gPad->SetLogx();

  TH1D *hchi = 0x0;
  const TString hns[]={"GausML", "CauchyML", "GausLS", "CauchyLS"};
  for(Int_t ii=0; ii<4; ii++){
    TH1D * htmp=(TH1D*)gDirectory->Get(Form("h2_%schi", hns[ii].Data()));
    if(htmp){
      if(hchi){
        printf("\nbad!! hchi already exist! %s\n", hchi->GetName()); exit(1);
      }
      else{
        hchi=htmp;
      }
    }
  }
  if(!hchi){
    printf("\nbad!! null hchi!!\n"); exit(1);
  }

  hchi->SetMaximum(7.5);
  hchi->SetMinimum(0);
  hchi->Draw();

  TF1 * f1=new TF1("f1","1",0,10);
  style::ResetStyle(f1);
  f1->SetLineColor(kRed);
  f1->SetLineWidth(2);
  f1->Draw("same");

  c0->Print(Form("outplot/momDescription%schi.png", outn.Data()));
}

void subFit(TTree * tt,  TH2D *hh, const TString yVar, const TString xVar, const TString cut, const Int_t varopt, const Double_t fitmin, const Double_t fitmax)
{
  const TString tag(hh->GetName());

  TCanvas *cc= new TCanvas;
  style::PadSetup(cc);

  const Double_t lm=0.15;
  Double_t rm=0.14;
  const Double_t tm=0.13;
  const Double_t bm=0.17;

  cc->SetLeftMargin(lm);
  cc->SetRightMargin(rm);
  cc->SetTopMargin(tm);
  cc->SetBottomMargin(bm);

  gStyle->SetOptTitle(1);
  gStyle->SetOptFit(1);
  gStyle->SetTitleW(0.95);
  gStyle->SetTitleX(0.3);
  gStyle->SetStatW(0.15);
  gStyle->SetStatX(1);

  //===========================

  TList * lout=new TList;
  lout->Add(hh);

  //======
  const TString plvar(Form("%s:%s>>%s",yVar.Data(), xVar.Data(), hh->GetName()));

  printf("\n\ncheck plvar %s || cut { %s }\n", plvar.Data(), cut.Data());
  tt->Draw(plvar,cut);

  hh->SetTitle(Form(";%s;%s",getTitleFromVar(xVar).Data(), getTitleFromVar(yVar).Data()));
  hh->Draw("colz");
  cc->Print(Form("outplot/%s2d.png", tag.Data()));

  //============================================================
  rm = 0.03;
  cc->SetRightMargin(rm);

  const Double_t thres = 20;
  TH1D *hrmsnor=0x0, *hrmsmpv=0x0, *hrmswid=0x0, *hrmsres=0x0, *hrmschi=0x0;
  NeutrinoTools::FitSlicesY(hh, hrmsnor, hrmsmpv, hrmswid, hrmsres, hrmschi, "RMS", thres, lout);

  TH1D *hnor=0x0, *hmpv=0x0, *hwid=0x0, *hres=0x0, *hchi=0x0;
  NeutrinoTools::FitSlicesY(hh, hnor, hmpv, hwid, hres, hchi, "CauchyML", thres, lout);

  hwid->SetMarkerStyle(20);
  hwid->SetMarkerSize(1);
  hwid->SetMarkerColor(kBlack);
  hwid->SetLineWidth(1);
  hwid->SetLineColor(kBlack);

  hrmswid->SetMarkerStyle(22);
  hrmswid->SetMarkerSize(1.5);
  hrmswid->SetMarkerColor(kRed);
  hrmswid->SetLineColor(kRed);
  hrmswid->SetLineWidth(1);

  TH1D * hfit = 0x0;
  TH1D * hrms = 0x0;
  if(varopt==0){
    hfit=hrmsmpv;
    hfit->SetMaximum(1.2);
  }
  else{
    hfit=hwid;
    hfit->SetMaximum(0.6);
    //hrms=hrmswid;
    //hrms->SetMaximum(0.6);
  }
  
  const Double_t hmax = hrms? hrms->GetBinContent(hrms->GetMaximumBin()) : hfit->GetBinContent(hfit->GetMaximumBin());
  //hfit->SetMaximum(hmax*1.3);
  hfit->SetMinimum(0);
  TF1 * ff=new TF1("ff","[0]/pow(x,[1])", fitmin, fitmax);
  ff->SetParameters(0.2,1);
  if(hrms){
    hrms->Fit(ff, "R");
  }
  else{
    hfit->Fit(ff, "R");
  }
  hfit->SetTitle(Form("%s  %s",hh->GetName(), ff->GetTitle()));
  lout->Add(ff);

  style::ResetStyle(hfit);
 

  if(hrms){
    style::ResetStyle(hrms);
    hrms->Draw();
  }

  hfit->Draw(hrms?"same":"");

  cc->Print(Form("outplot/%svaropt%d.png", tag.Data(), varopt));

  //========

  TFile * fout=new TFile(Form("outplot/fitslice%s.root", tag.Data()),"recreate");
  lout->Write();
  fout->Save();
  fout->Close();
  delete fout;
}

void fitLeverArm(const TString fin, const Int_t mode)
{
  TFile::Open(fin);
  TTree *tt = tree;

  gROOT->LoadMacro("style.cxx+");

  style::SetGlobalStyle();

  style::fgkYTitleOffset = 1.4;
  style::fgkXTitleOffset = 1.4;

  //===================================================
  const TString basecut = "fMuonNRec>0 && fMuonNSim>0 && fMuonRec->P()<3.5";

  if(mode < 10){
    const Double_t nhitmin=55+mode*5;
    const Double_t nhitmax = nhitmin+5;
    const TString tag = Form("nhit%.0f%.0f", nhitmin, nhitmax);
    const TString nhitcut = basecut+Form(" &&  fMuonRecNhits[0]> %.0f && fMuonRecNhits[0]<%.0f", nhitmin, nhitmax);
    //subFit(tt, new TH2D(tag+"dd","", 12,0, 4, 60, -2, 2),  "(1/fMuonRec->P()-1/fMuonSim->P())/(1/fMuonRec->Pt()-1/fMuonSim->Pt())",   "fMuonRecFlightPath->Pt()", nhitcut, 1, 0.7, 3.1);
    subFit(tt, new TH2D(tag+"PToverP","", 12,0, 4, 110, 0, 1.1), "fMuonSim->Pt()/fMuonSim->P()",      "fMuonRecFlightPath->Pt()", nhitcut, 0, 0.7, 3.1);
    subFit(tt, new TH2D(tag+"dpT","", 12,0, 4, 60, -1.2, 1.2),   "1/fMuonRec->Pt()-1/fMuonSim->Pt()", "fMuonRecFlightPath->Pt()", nhitcut, 1, 0.7, 3.1);
    subFit(tt, new TH2D(tag+"dMom","", 12,0, 4, 60, -1.2, 1.2),  "1/fMuonRec->P()-1/fMuonSim->P()",   "fMuonRecFlightPath->Pt()", nhitcut, 1, 0.7, 3.1);
  }
  else{
    mode -=10;
    const Double_t larmmin=0.5+mode*0.5;
    const Double_t larmmax = larmmin+0.5;
    const TString tag = Form("larm%03.0f%03.0f", larmmin*100, larmmax*100);
    const TString larmcut = basecut+Form(" &&  fMuonRecFlightPath->Pt()> %.2f && fMuonRecFlightPath->Pt()<%.2f", larmmin, larmmax);
    
    subFit(tt, new TH2D(tag+"PToverP","", 10, 0, 100, 110, 0, 1.1), "fMuonSim->Pt()/fMuonSim->P()",      "fMuonRecNhits[0]", larmcut, 0, 30, 80);
    subFit(tt, new TH2D(tag+"dpT","", 10, 0, 100, 60, -1.2, 1.2),   "1/fMuonRec->Pt()-1/fMuonSim->Pt()", "fMuonRecNhits[0]", larmcut, 1, 30, 80);
    subFit(tt, new TH2D(tag+"dMom","", 10, 0, 100, 60, -1.2, 1.2),  "1/fMuonRec->P()-1/fMuonSim->P()",   "fMuonRecNhits[0]", larmcut, 1, 30, 80);
  }
}

cmpMomMuon()
{
  gStyle->SetOptStat(0);

  c1=new TCanvas;

  gPad->SetLogx();
  gPad->SetLogz();

  hh=new TH2D("hh","",100,0.05,2,300,0,3);NeutrinoTools::BinLog(hh->GetXaxis());
  
  tree->Draw("fTmpMom/fMuonRec->P():fMuonRecPt->Mag()>>hh","fMuonNRec>0 && ( fMuonRec->P()<3.5 && fMuonRecNhits[0]>60 && fMuonRecFlightPath->Pt()>1.5  ) && ( fMuonRecNhits[0]<80 && fMuonRecFlightPath->Mag()>3.0 )","colz"); 
  h2=NeutrinoTools::NormalHist(hh,2,1); h2->Draw("colz");
  h2->SetTitle("mom'/mom:pT, mom'=reconstructed with muon hypo. after QA cuts");
  h2->GetYaxis()->SetRangeUser(0.5,1.5);
  h2->Draw("colz");
  c1->Print("outplot/momMuon1.png");

  c1=new TCanvas;
  gPad->SetLogx();
  gPad->SetLogz();

  tree->Draw("fTmpMom/fMuonRec->P():fMuonRecPt->Mag()>>hh","fMuonNRec>0","colz"); 
  h2=NeutrinoTools::NormalHist(hh,2,1); h2->Draw("colz");
  h2->SetTitle("mom'/mom:pT, mom'=reconstructed with muon hypo. no QA cuts");
  h2->GetYaxis()->SetRangeUser(0.5,1.5);
  h2->Draw("colz");
  c1->Print("outplot/momMuon2.png");
}
